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Abstract 


Temporary capture efficiency is studied in the framework of the circular restricted three-body problem in two steps. 
First, a non-uniform distribution of test particles around the secondary’s orbit is obtained by fully accounting the 
secondary’s gravitational influence. Second, the capture efficiency is computed based on the non-uniform 
distribution. Several factors influencing the result are discussed. By studying the capture efficiency in the circular 
restricted three-body problem of different mass ratios, a power-law relation between the capture efficiency (p) and 


the mass ratio (u) is established, which is given by px0.27 x wu 


0-53 within the range of 3.0035 x 


10 < u < 3.0034 x 10%. Taking the Sun—Earth system as an example, the influence from the orbit 
eccentricity of the secondary on the non-uniform distribution and the capture efficiency is studied. Our studies 
find that the secondary’s orbit eccentricity has a negative influence on the capture efficiency. 
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1. Introduction 


Gravitational capture is an important dynamical phenomenon in 
the evolution of the solar system, which has engaged the curiosity 
of many researchers. It usually occurs when a celestial body 
approaches a massive planet at a relatively low velocity during a 
close encounter, leading to a transition from a heliocentric orbit to 
a planetocentric one (Heppenheimer & Porco 1977). During this 
process, the orbit may be significantly perturbed from a Keplerian 
orbit about the Sun as the gravitational perturbation from the close 
encounter planet becomes dominant. 

Hopf (1930) and Tanikawa (1983) demonstrate that the 
probability of permanent capture in a purely gravitational 
environment is zero, yet temporary capture is always possible. 
When combined with some dissipative mechanisms, the tempor- 
ary capture can become permanent. The dissipative mechanisms 
can be gas drag (e.g., Hunten 1979; Pollack et al. 1979; Cordeiro 
et al. 1999; Vieira Neto et al. 2009), a gradual increase in the mass 
of the central planet (e.g., Heppenheimer & Porco 1977; Vieira 
Neto et al. 2004), and collisions (e.g., Colombo & Franklin 1971). 
The capture can also become permanent when a binary asteroid 
system encounters a planet, in which case one of the bodies is 
permanently attached to the planet and the other escapes due to 
the exchange of orbital energies (e.g., Agnor & Hamilton 2006). 
Due to these mechanisms, the capture process provides 
researchers with an opportunity to explain the origin and evolution 
of irregular satellites around four outer solar system planets (e.g., 
Nesvorny et al. 2007, 2014). 

The capture process also arouses interest from the aerospace 
community as it offers potential fuel savings for interplanetary 


spacecraft by strategically utilizing the capture mechanism to 
achieve a nominal orbit around a planet. Belbruno (1987) first 
proposed the concept of Weak Stability Boundary and 
successfully applied it to the Japanese Hiten spacecraft 
(Belbruno & Miller 1993). Topputo & Belbruno (2009), Luo 
et al. (2014), Luo & Topputo (2017, 2021), Luo (2020) studied 
the gravity-assisted capture in different systems. Studies have 
also been carried out on the application of the capture 
mechanism to retrieve asteroids into the near-Earth space 
(e.g., Baoyin et al. 2010). 

Solving the capture problem analytically is challenging. The 
capture occurs when celestial bodies encounter the planet at a 
low relative speed. In such a case, the dynamics are usually 
chaotic. Murison (1989) and Cordeiro et al. (1999) pointed out 
the complexity of the capture process in the circular restricted 
three-body problem (hereafter referred to as CRTBP), since the 
capture structure in the phase space is partially fractal-like and 
self-similar. Researchers therefore generally prefer to use 
numerical methods to study the capture process, and usually 
make some reasonable assumptions to simplify the situation 
and reduce the computational difficulty. Carusi & Pozzi (1978), 
Carusi & Valsecchi (1980) worked on the capture phenomenon 
of objects located near Jupiter by setting the initial velocities of 
particles tangential to Jupiter’s orbit. Huang & Innanen (1983) 
studied the stability of retrograde jovicentric orbits to find the 
relationship between the length of the capture time and the 
periodic orbits using Henon’s diagram (Henon 1970). The 
mirror theorem (Roy & Ovenden 1955) is also used to simplify 
the problem whose nature is time reversible and symmetric in 
the CRTBP. In the study of de Almeida Prado & Vieira Neto 
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(2006), a polar mesh grid of distance and angle around the 
capturing planet was adopted. They used the reversibility of 
time in the CRTBP and integrated the particles backwards until 
they left the vicinity of the planet. Brunini (1996) explored the 
capture conditions for particles by placing particles on a fixed 
line perpendicular to the x-axis. He changed the angle a 
between the velocity and the x-axis and computed the capture 
domain in the phase space of (y, a). In addition to these studies, 
some others also focus on the length of the temporary capture 
time (e.g., Vieira Neto & Winter 2001; Winter & Vieira 
Neto 2001; Fedorets et al. 2017). The capture process is also an 
important part of the transport of matter in interplanetary space. 
It is the basis for calculating the flux of matter orbiting and 
impacting the planet. 

A practical concern of this study is the computation of the 
temporary capture efficiency within a general CRTBP system, 
taking into account the gravitational influence of the secondary 
(hereafter referred to as P2 and the primary referred to as P1). It 
represents the probability of a massless particle being 
temporarily captured by P2 within a prescribed duration. In 
our work, the capture efficiency serves as a metric for assessing 
the capability of P2 to capture nearby particles in the CRTBP. 
The basic method of computing the capture efficiency in our 
work is the one employed by Granvik et al. (2012) which 
computes the temporary capture efficiency of near-Earth 
objects (hereafter referred to as NEOs). It involves the 
following steps: 


1. Choose a ring region enveloping the Earth’s orbit. The 
ring is described by the semimajor axis, orbit eccentricity 
and inclination, and is discretized by uniformly distrib- 
uted bins inside it. 

2. Inside each bin of the above ring region, generate N; 
random test particles with random values of longitude of 
the ascending node Q, argument of periapsis w, and mean 
anomaly M in the range of [0, 27]; 

3. Choose a circumsphere centered at the Earth inside the 
ring region, and count the number of test particles inside 
the circumsphere (denoted as N>); 

4. Integrate orbits of the N3 particles initially inside the 
circumsphere within a prescribed time T and count the 
number of particles that are temporary captured by the 
Earth (denoted as N3); 

5. The capture efficiency of the bin is computed as p(a, e, i) = 
N3/N1. 


The total capture efficiency of near-Earth objects (NEOs) by 
the Earth is then obtained by multiplying the orbital element 
distribution R(a, e, i) of NEOs (Bottke et al. 2002) with the 
capture efficiency in each bin and then summing over all 
the bins. 

In the current study, we use the same method to compute the 
capture efficiency, but focus on the general CRTBP model. The 
orbital element distribution function R(a, e, i) is generated by 
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fully considering the resonance effects between P2 and the 
particles. Moreover, we made following contributions: (1) A 
backward integration method is proposed to identify the ring 
region enveloping the orbit of P2. (2) Influence of the 
circumsphere size on the capture efficiency is investigated, 
and appropriate values are recommended. (3) By computing the 
capture efficiency of CRTBPs with different mass ratios, a 
power-law relation between the capture efficiency and the mass 
parameter of CRTBP is established. (4) Influence of P2’s orbit 
eccentricity on the capture efficiency is studied. 

The paper is organized as follows. Section 2 is devoted to the 
methods used in the study, including the way to obtain the non- 
uniform steady-state distribution which is influenced by the 
resonances with P2, the way to choose the ring region, the 
circumsphere size, and the integration time, and the way to 
compute the capture efficiency. Section 3 presents the main 
results of our numerical experiments. A power-law relation 
between the capture efficiency p and the mass ratio ju is given. 
Section 4 discusses the influence of P2’s orbit eccentricity on 
the capture efficiency. Section 5 concludes the study. 


2. Methodology 
2.1. Capture Efficiency Computation 


The method to compute the capture efficiency is the one used 
by Granvik et al. (2012), which mainly contains the following 
steps: 


1. The first step is to choose an initial sampling region 
enveloping P2’s orbit. The sampling region is character- 
ized by [@min» Amax] x [emin; max] x Limin» imax]. It is 
discretized into a series of bins with uniform step sizes 
along the semimajor axis, orbit eccentricity, and orbit 
inclination. 

2. Inside each bin, generate a number of test particles whose 
Q, w, and M are randomly distributed within the range of 
[0, 27]. The total number of test particles in each bin is 
denoted as Mgen. 

3. Choose a circumsphere inside the ring region, with its 
center at P2. Select test particles that are inside the 
circumsphere. Denote the number of the selected test 
particles as Méiter- 

4. Integrate orbits of the selected test particles and count the 
number of test particles that can be temporarily captured 
within an integration time 67. Denote the number of 
captured particles as M¢ap. 

5. The capture efficiency of this bin is computed by the 
following formula: 


Mcap Myitter 


pla, e, i) = l (1) 
Miter Meen 


6. Suppose that the distribution of the test particles around 
P2 is in a steady-state, and is described by a function 
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R(a, e, i). The total capture efficiency is computed by the 
following formula: 


Psum = Nr. e, i)R(a, É, i)dadedi 
Noin 
= >` pla, e, i)R(a, e, i) (2) 


i=l 
in which Npjn is the total number of discretized bins. 


In their study, Granvik et al. (2012) used the steady-state 
distribution of NEOs Rygo(a, e, i) given by Bottke et al. 
(2002). They also give recommendations of the ring size, 
circumsphere size and integration time for the Sun—Earth 
system. However, for the purpose of the current study which 
focuses on the general CRTBP model, we have to find these 
settings by ourselves before applying the above method. In the 
following subsections, we explain more details. 


2.2. Capture Criterion and Boundary 


The temporary capture criterion is rather a controversial but 
inevitable issue as there is no precise dynamical constraint on 
whether a celestial body is captured. Everhart (1973) judged 
whether a celestial body is captured by monitoring if the 
osculating planetocentric orbital elements become elliptical. 
However, Fedorets et al. (2017) mentioned that quasi-satellites 
may have an osculating orbit eccentricity less than 1 with 
respect to the planet, while they are not dynamically bound to 
the planet. Carusi & Pozzi (1978) and Gladman et al. (1995) 
checked the distance of the particles from the planet. 
Yamakawa et al. (1992) related the capture efficiency with 
the two-body energy and assumed that the capture occurs when 
the energy of the celestial body relative to the planet becomes 
negative. In Granvik et al. (2012), the same definition as the 
one used in Kary & Dones (1996) is adopted which includes 
two criteria: 


1. The distance of the celestial body away from P2 is less 
than three Hill Radii of P2. 

2. The two-body Keplerian energy of the celestial body with 
respect to P2 is less than 0. 


Once the particles in the CRTBP system simultaneously satisfy 
the above two conditions, it can be considered to be 
temporarily captured by P2. In this study, we take the same 
two conditions. In the discussion section, we will discuss why 
we persist using the value of three Hill Radii. 


2.3. Initial Sampling Region 


One issue to be addressed is how to choose the initial 
sampling region for a given CRTBP system. Generally, the size 
of the sampling region should increase with the mass 
parameter, as the gravitational influence region of P2 expands. 
We need to specify the initial sampling region in which the 
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Capture efficiency analysis in the CRTBP 


Initial sampling region 


Figure 1. The schematic diagram of the initial sampling region and the 
circumsphere. The radius Ffiter corresponds to the radius of the circumsphere 
and should be inside the initial sampling region. The radius rin and Fmax 
represent the inner and the outer boundary of the ring-shaped initial sampling 
region, which approximately equal the minimum periapsis distance and the 
maximum apoapsis distance of the particles. 


Figure 2. The schematic diagram of the particles being captured and passing 
the circle with a radius of three Hill Radii of P2. Particles are assumed to have 
zero energy when they pass through this circle. v is the inbound velocity of the 
particle at the circle. 


gravity of P2 has an obvious influence and the particles have a 
good chance of being captured (see Figure 1). In this 
subsection, we propose a backward integration method to 
identify the initial sampling region. 

A necessary condition for temporary capture is that particles 
have a low relative velocity at close encounters with P2. This 
implies that the particles are on orbits characterized by low 
eccentricities and inclinations with respect to P2 (Qi & 
Qiao 2023). As a result, we can use the planar CRTBP model 
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Figure 3. The range of the initial sampling region determined by the backward integration method is shown. The mass parameter is the one of the Sun—Earth system 
(ue 3.0035 x 1076). The left frame shows the histogram of the reachable bin in the (a, e) space. The right frame shows the identified initial sampling region (the 


smaller square) and the one used by Granvik et al. (2012) (the larger square). 


to identify the initial sampling region in the (a, e) space. An 
intuitive idea is that in the planar CRTBP, trajectories of the 
captured particles should pass through a circle surrounding P2 
with an instantaneous incoming velocity v which forms an 
angle 0 > 7/2 with respect to its position vector, as depicted in 
Figure 2. By generating points on the circle and performing a 
backward integration of them from the “capture boundary” 
(three Hill Radii which is the same as that used in Granvik et al. 
2012), we can obtain the region in the space of (a, e) where the 
particles originate. Since it is possible for the particles to enter 
the three Hill Radii region at any point on the sphere, the entire 
circle needs to be surveyed. The energy of the particle with 
respect to P2 is set to zero as it is also one of the temporary 
capture criteria used by us. Hence the value of |v| can be solved 
by this condition. 

Using the mass parameter of the Sun—Earth system as an 
example (ug 7% 3.0035 x 1076), 100 points are selected uni- 
formly at a distance of three Hill Radii from P2. At each point, 
we create 100 random directions for the inbound velocity 
vector v ensuring that the constraint @>7/2 is met. 
Consequently, there are a total of 100 x 100 orbits. These 
orbits are backwards integrated for 100 units of dimensionless 
time in the CRTBP. For each orbit, we record its minimum and 
maximum orbit semimajor axes, as well as its maximum orbit 
eccentricity. The left frame of Figure 3 shows the statistical 
result. The plane spanned by the semimajor axis and the orbit 
eccentricity is divided into a series of bins (the resolution is 
Aa = 0.01 au, Ae = 0.01), and the relative weight of each bin 
is computed by counting the number of particles in it. If we 
neglect the bins that contain particles less than 10% of the 
maximum one, we get the initial sampling region (denoted by 
the smaller square) in the (a, e) space presented in the right 
frame of Figure 3. The larger square represents the initial 
sampling region used in Granvik et al. (2012). The two regions 


Table 1 

Parameters Used in the Numerical Experiments 
Mass Para- Qmin Dias 
meter (4g)? (au)? (au)? ema? Noin? Men Mos 
1 0.87 1.15 0.06 840 8.4 x 10° 8.4 x 10* 
2 0.85 1.20 0.08 1400 1.4 x 107 1.4 x 10° 
3 0.83 1.23 0.09 1800 1.8 x 107 1.8 x 10° 
4 0.81 1.26 0.10 2250 225x110’ 2.25 x 10° 
5 0.80 1.29 O11 2695 2.695 x10” 2.695 x 10° 
6 0.79 131 0.12 3120 312x110’ 3.12x 10° 
7 0.78 1.33 0.12 3300 3.3 x 107 3.3 x 10° 
8 0.77 135 013 3770 3.77x10’ 3.77 x 10° 
9 0.76 1.37 0.14 4270 427x10 4.27x 10° 
10 0.76 1.38 0.14 4340 4.34 x 107)  4.34x 10° 
Notes. 


è The mass parameter of the CRTBP systems (ue ~ 3.0035 x 107°). 
b en : ; 
The range of the initial sampling region. 
© The total number of the discretized bins. 
d The total number of the particles for computing the capture efficiency of each 
bin. 
° The total number of the particles for computing the non-uniform steady-state 
orbit distribution. 


agree with each other quite well in the semimajor axis, but not 
in the orbit eccentricity. The upper boundary of Granvik et al. 
(2012) is about twice the upper boundary of ours. 

With the increase of P2’s mass, the size of the initial 
sampling region also increases. Figure 4 shows the range of 
initial sampling region when the mass parameter of the system 
is t=10pg, The range apparently enlarges compared to that in 
Figure 3. Table 1 lists the relevant parameters used in our work. 
The bin resolution is set to Aa=0.01 au, Ae=0.01. As we 
have mentioned, the capture process mainly happens for planar 
orbits, so in our study the orbit inclination is restricted as 
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Figure 4. Same as Figure 3 but for a CRTBP with mass parameter of u = 10pUp. 


i € [0°, 2°5] which is in accordance with Granvik et al. (2012), 
and the bin size in the orbit inclination direction is 
Ai=0.5°.=0°5. We generate 10° particles in each bin to 
compute the capture efficiency p(a, e, i) of this bin. 


2.4. Circumsphere Size and Integration Time 


It is unnecessary to integrate all the test particles in the initial 
sampling region. Only test particles close to P2 have a chance 
to be temporarily captured within a prescribed integration time. 
Numerically integrating all the test particles is time-inefficient 
and unnecessary. A clever way to avoid integrating all the test 
particles is to choose a circumsphere around P2 to filter those 
that are close to P2 as integration candidates in Granvik et al. 
(2012). We use the same method in our work. The schematic 
diagram of the circumsphere is shown in Figure 1. These 
selected candidates are further integrated for a prescribed time 
to judge whether they are captured or not. The radius of the 
circumsphere is denoted as /fiter and the length of the 
integration time is denoted as 67. The final computed capture 
efficiency p is therefore a function dependent on these two 
parameters. It is only an approximation of the accurate capture 
efficiency which should be obtained by integrating all the test 
particles in the initial sampling region for a long enough 
integration time. 

A natural expectation is that as the circumsphere size 
increases, the capture efficiency gradually approaches its 
accurate value. We discretize the circumsphere size in units 
of Hill Radius of P2 as Ffiter = k - Ry(k € N). The following 
relationship should hold: 


om p(k + Ry, 6T) = p(6T). (3) 
For the purpose of practical computation, we need to select an 


appropriate value of k. One remark is that the circumsphere 
should be within the initial sampling region. This puts an upper 


limit for the value of k, denoted as kmax. For each k € [0, kmaxl, 
we compute the capture efficiency p(k). The size selection 
should be large enough in order not to miss the particles that 
can be captured during a specific integration time. Starting from 
smaller values, we select the k value once the following 
criterion is met: 


p(k: Ru, ôT) — p(k — 1) - Ru, ôT) 


Ap(k) = 
pg p(k — 1) - Rw, 67) 


< 10% (4) 
where Ap(k) stands for the relative increment in capture 
efficiency for a specific circumsphere size compared to the 
former one. If the criterion cannot be satisfied unless k = kmax, 
we Select kmax aS our k value. 

As for the integration time 67, it is obvious that the longer it 
is, the better the computed result is. However, considering the 
practical computation cost, we need to put an upper limit to it. 
In the study of Granvik et al. (2012), the integration time is 
chosen as 2000 days. In our work, we choose 2 yr for the Sun— 
Earth system. We will show that the option of this integration 
time is appropriate. 


2.5. Steady-state Orbit Distribution 


According to Equation (2), to get the actual capture 
efficiency, we need a distribution of test particles around 
P2’s orbit. In the work of Granvik et al. (2012), an NEO orbit 
distribution model Rygo(ap, €n, in) from Bottke et al. (2002) is 
utilized. For the general analysis in this study, we do not have 
the information and we have to generate it by ourselves. 

In reality, the distribution of small bodies around a planet’s 
orbit is the result of complicated evolution processes during 
which the source regions, the gravitational effects, and even 
some non-gravitational effects (such as collisions and thermal 
effects) may have their contributions. For the simple CRTBP 
model considered in this study, P2’s gravity is the main 
contribution for a non-uniform steady-state distribution. It is 
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Figure 5. Left frame: One example bin in the two-dimensional phase space spanned by a and e. Right frame: The number of particles in one example bin (ID number 


is 36) as a function of the evolutionary epoch jj. 


obvious that due to orbital resonances between P2’s orbit and 
the particles’ orbit, the steady-state distribution around P2’s 
orbit is non-uniform. To get the non-uniform orbital element 
distribution, we start with an initial uniform distribution in the 
initial sampling region. Take the mass parameter of the Sun- 
Earth system as an example. The range of this region is 
a€ [0.87, 1.15] au, e € [0, 0.06], ¿€ [0°, 2°5], and the bin 
resolution is Aa = 0.01 au, Ae = 0.01, Ai=0°5. The angular 
elements (Q, w, M) are uniformly distributed in the range of 
[0, 27]. 

A time-slicing method is used. The nature of the method is to 
simulate a dynamically steady-state population by means of 
time discretization. We assume that the orbital element 
distribution of the particles around P2’s orbit is already in a 
steady-state. That is to say, influenced by P2’s gravity, at any 
given time, the number of particles leaving the sampling region 
is equal to the number of particles entering this particular 
region. At any given observation epoch T and in a specific 
region of the orbital element space spanned by (a, e, i), we 
observe the population of the particles. The time is discretized 
from the initial epoch (fọ =0) to the epoch T into N equal 
intervals, so that we have a series of time epochs as t; =i- AT 
(i=0--- N) where the time interval AT = T/N. The particles 
observed are a mixture of some old particles that were in this 
region before epoch T, and some new particles that have just 
entered this particular region at epoch ty = T. 

To simulate this process, we perform the following 
numerical simulations: 


1. Since we are dealing with a continuous system, we first 
discretize the time into a series of epochs t;=i- AT. The 


epoch T is the total evolution time and also the 
observation epoch at which we “observe” the orbital 
element distribution, and N is the number of time 
intervals in the total evolution process, i.e., N= T/AT. 
We then divide the phase space of the orbital elements 
into a set of bins according to certain resolution. Using 
the two-dimensional phase space spanned by a and e as 
an example, the left frame of Figure 5 shows one example 
of the binned orbital element space. 


. At observation epoch T, we count the number of particles 


in each bin. As we have mentioned, the particles in each 
bin include those that have only entered the bin at the 
epoch T, and those that have entered the bin at an earlier 
epoch ¢; < T and still remain in the bin after an evolution 
interval of T — t;. Denote the number of these particles as 
MË ) where the subscript i is the group number of 
particles and the subscript j is the ID of the binned orbital 
element space. The total number of the particles in the bin 
(j) is therefore MY? = DA 0 MY, which is the sum of all 
previous particles that are in the bin. 


. Suppose that the particles in all bins have an equal chance 


of evolution. To approximate this equal chance, we give 
each bin the same number of initial test particles to 
simulate the initial uniform distribution and track their 
evolution over time. Starting from the initial epoch tọ = 0, 
we count the total number of particles from all bins that 
can enter the bin(j) at each epoch t; Of course, the 
number changes with the epoch t; It is necessary to 
distinguish between the epoch ft; which stands for the 
group number of particles starting the evolution and the 
integration epoch which is the evolution time in this 
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context. If we integrate the test particles for a total time of 
T and denote the integration time as f, the number at f; 
corresponds to the particles that have entered the bin at 
epoch ft,_; and remain in the bin after an interval of 
T — ty_;. One example bin is shown in the right frame of 
Figure 5. According to the curve, it is reasonable to say 
that the value MẸ’, at the evolution epoch j is the 
number of particles that are in the bin(/) after a time 
interval of T — ty_;. For example, when i = 0, the number 
M G ) at epoch fy indicates the initial test particles that have 
not started their evolution yet, while if i= N, the number 
M$” at epoch iy indicates the final test particles in the bin 
after an evolution time of T. 


The orbital element distribution at the epoch T is thus 


Ma, e, i) 


R(a, e, i, T) = ue 


, M*=S°M% a, e,i) (5) 
where M“ is the total number of the particles at “observation” 
epoch T in bin(j), and MË% is the total number of the particles 
at “observation” epoch T in all bins. Since a steady-state orbit 
distribution is assumed, the above process should be indepen- 
dent of the observation epoch T (also the evolution integration 
time) and R(a, e, i) = R(a, e, i, T). Intuitively, the longer the 
integration time T, the larger the total number of the test 
particles for the evolution M$, the smaller the time interval 
AT, and the higher the resolution of the bins, the closer the 
resulting orbital distribution will be to the true case. 

In our study, the evolution integration time is set as 
T=10° yr, and the time interval for “observing” the orbit 
distribution is 104 yr. The bin interval is mentioned above as 
Aa=0.01 au, Ae=0.01, Ai=0°5. To carry out the above 
numerical integration process, we generate 100 particles in 
each bin to track their orbit evolution. Denote the total number 
of particles generated to obtain the non-uniform steady-state 
distribution as M$. According to Table 1, within 840 bins, a 
total of 8.4 x 10* particles are generated. It is necessary to 
distinguish the number of Mow and Mn as they correspond to 
two different numerical processes. 

One remark is that the assumption that particles in all bins 
have an equal chance of evolution is inappropriate. Given that 
the orbit distribution near P2 is supposed to be in a steady-state 
and non-uniform, there is no doubt that fewer particles should 
be in a bin if it has a lower weight R’, and vice versa. To take 
this non-uniformity into account, a basic idea is to repeat the 
above process after obtaining the approximation of the orbital 
element distribution R(a, e, i) based on the initial uniform 
distribution. This time, the number of initial particles in each 
bin for a new evolution process should be chosen according to 
R(a, e, i), which is also the weight of the bin. However, we did 
not iterate in this study. In reality, a steady-state distribution is 
a balance between external sources which keep sending 
particles into the initial sampling region and the particles that 


Miao & Hou 


are already in the initial sampling region. Under the sole 
influence of P2’s gravity, the region with a high weight means 
that orbits of particles in these regions are more stable. On the 
other hand, particles from external sources are harder to get into 
these regions due to their stability. As a result, when 
considering external sources, these regions do not necessarily 
have a higher weight. Since in our general analysis, we do not 
have information of external sources, so we prefer not to iterate 
because otherwise stable regions will gain a higher and higher 
weight after each iteration. We emphasize that since we only 
consider P2’s gravity in our work, the capture efficiency only 
reflects the ability of P2 to temporarily capture particles, but is 
not the real capture efficiency in real physical systems in which 
real steady-state distribution should be used. 


2.6. Numerical Integrators 


In our work, two integration algorithms are employed. To 
compute the non-uniform orbital element distribution intro- 
duced in Section 2.5, a hybrid symplectic integrator package, 
Mercury6 written by Chambers (1999) is used to compute the 
dynamical evolution of particles. The hybrid symplectic 
algorithm can accurately handle close encounters by predicting 
and using a conventional Bulirsch—Stoer extrapolation inte- 
grator to step over it. To compute the capture efficiency, a 
Runge—Kutta—Felberg 7(8) integrator is used to follow the 
orbits of the particles within the circumsphere. A particle is 
removed once it impacts P2 or P1 if its instantaneous P2-centric 
distance is less than the radius of P2 or its Pl-centric distance is 
less than the radius of P1. The radius of P1 is always assumed 
to be the same as the Sun, while the radius of P2 increases with 
mass ratio by keeping the density the same as that of the Earth. 


3. Results 


In this section, we first take the mass parameter of the Sun— 
Earth system (upg 3.0035 x 1076) as an example to display 
the results. After that, the relationship between the capture 
efficiency p and the mass parameter jz is shown. 


3.1. Capture Efficiency of Each Bin 


According to Table 1, in each bin 104 particles are uniformly 
generated within the initial sampling region. The range of 
initial sampling region is: a € [0.87, 1.15] au, e € [0, 0.06], 
i€[0°, 2°5] and the bin resolution is: Aa=0.01 au, 
Ae = 0.01, Ai=0°5. 

As we have mentioned in Section 2.4, it is necessary to first 
check the influence of the filtering circumsphere size on the 
capture efficiency, and choose a proper value for circumsphere 
size and integration time. Figure 6 displays the capture 
efficiency for different circumsphere sizes after the integration, 
as well as the relative increment in the capture efficiency 
between adjacent two circumsphere sizes defined by 
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Figure 6. The red line stands for the capture efficiency at each circumsphere 
size based on the uniform distribution. The yellow line stands for the relative 
increment in the capture efficiency. The black dashed line is the 10% criterion. 
At about 7 Hill radii, the relative increment is smaller than 10%. 


Equation (4). It seems that when k= 7 Hill radii, the relative 
increment of the capture efficiency is less than 10%. This is 
larger than the size of the circumsphere taken by Granvik et al. 
(2012). One remark is that if we set the truncation criterion as 
20%, the value of k is around 4-5, which is the value taken by 
Granvik et al. (2012). 

Fixing the circumsphere as seven Hill Radii, Figure 7 shows 
the relationship between the growth of capture efficiency and 
the integration time. An obvious trend is that the capture 
efficiency gradually approaches a fixed critical value when the 
integration time is long enough. It can be seen that when the 
particles are integrated for 2 yr, the growth of capture efficiency 
almost ceases. Therefore, a prescribed 2 yr integration time in 
the Sun—Earth system is long enough. 

Based on the chosen initial sampling region, the circum- 
sphere size, and the integration time, the capture efficiency of 
each bin can be computed. According to Equation (1), the map 
of capture efficiency p(a, e, i) in each bin is illustrated in the 
left frame of Figure 8. 


3.2. Steady-state Distribution 


To numerically simulate the long-term evolution of the 
whole population, a total number of 8.4 x 10* particles are 
integrated for T= 10° yr, as listed in Table 1. The time interval 
for observing the state of orbit distribution is AT= 10* yr. 
Figure 9 shows the final steady-state distribution of the test 
particles in the orbital element (a, e) space through the 
numerical method in Section 3.2. Due to the gravitational 
influence of P2, the distribution is non-uniform. The bins have 
a higher weight if the color is lighter. The cumulative ridge in 
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Figure 7. The blue line stands for the cumulative curve of the capture 
efficiency with respect to the integration time. The red line stands for the 
relative increment in the capture efficiency. The circumsphere size is set to 7 
Hill Radii. After 2 yr of integration, the growth in capture efficiency has almost 
ceased. A prescribed 2 yr integration time to compute the capture efficiency is 
long enough. 


the left frame of Figure 9 is located where the semimajor axis 
a= l au. Particles accumulate here because they are on orbits 
like P2 and are trapped in the 1: 1 resonance. The region with 
lower weight close to 1 au corresponds to P2-orbit crossers 
with low orbital eccentricities. Close encounters between P2 
and the particles are common here. P2 wipes out the particles 
here and causes a significant depletion in this region. Particles 
with semimajor axes further away from | au and small orbit 
eccentricities have little chance to encounter P2, so the bins in 
this region have a higher distribution weight. When the bin 
resolution is enhanced to Aa = 0.0005 au, Ae = 0.0005 (see 
the right frame of Figure 9), we can observe some fine 
structures at specific positions. We believe that these structures 
are closely related to orbital resonances with P2, which are very 
close to the 1: 1 one (Pan & Hou 2022; Tan et al. 2023). 

According to Equation (2), the capture efficiency p(a, e, i) of 
each bin is multiplied with the bin weight R(a, e, i). The result 
is illustrated in the right frame of Figure 8. Based on the non- 
uniform steady-state distribution and the capture efficiency of 
each bin, the capture efficiency is computed as ~2.930 x 10-4, 
which is only ~50.49% of the value if using a uniform steady- 
state distribution. 


3.3. Relationship between the Capture Efficiency and the 
Mass Ratio 


In this section, we study the relation between the capture 
efficiency p and the mass ratio u. The values of the mass 
parameters are listed in Table | and the corresponding capture 
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Figure 8. The left frame is the capture efficiency p(a, e, i) of each bin for a mass parameter of the Sun—Earth system. The right frame is the modified capture efficiency 
distribution based on the non-uniform steady-state distribution. The colorbar means the capture efficiency of the bin. 
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Figure 9. The steady-state distribution of the particles near P2’s orbit for the system with js ~ 3.0035 x 1076. The colorbar is the relative weight of the bin. (1) The 
left frame: The distribution in the (a, e) space with the bin resolution as Aa = 0.01 au, Ae = 0.01. The red line stands for test particles with apoapsis distance of 1 au 
and the black one stands for test particles with periapsis distance of 1 au; (2) The right frame: The distribution in the (a, e) space with the resolution as 
Aa = 0.0005 au, Ae = 0.0005. The red dashed lines stand for the positions of the orbital resonances with P2, which match the computed fine structures well. 


efficiency values are computed. The procedure of computing 
the capture efficiency for different mass ratios u is the same as 
that for the Sun—Earth system described in Sections 3.1 and 
3.2. The sizes of the initial sampling regions and the total 
numbers of test particles used in the numerical simulations are 
also listed in Table 1. The results of the capture efficiency for 
different mass ratios based on different orbit distributions are 
displayed in Figure 10. The red dots and fitted short-dash line 
correspond to the results of the uniform steady-state distribu- 
tion of test particles, i.e., P2’s gravitational influence on the 
orbital element distribution is ignored. The green dots and fitted 
dotted—dashed line correspond to the results of the non-uniform 
distribution of test particles, which is obtained by fully 
considering P2’s gravitational influence. Obviously, the capture 
efficiency of the non-uniform distribution is smaller than the 
result of assuming a uniform distribution. For the non-uniform 
steady-state which considers the gravitational influence of P2, a 
power-law relation between the capture efficiency p and the 
mass ratio p is fitted as 


p x 0.27 x p> 


in which 3.0035 x 1076 < u < 3.0034 x 107°. 


4. Discussion 
4.1. Capture Boundary 


In Section 2.2, if we refer to the P2-centric distance in 
Criterion 1 as “capture boundary”, it is natural to wonder 
whether the final computed result would be significantly 
changed when this value is altered. In previous work, Granvik 
et al. (2012) and Kary & Dones (1996) did not explain the 
reason why three Hill Radii is chosen as the boundary to judge 
whether a particle is captured by the planet. Intuitively, the 
farther the particles are away from the planet, the lower the 
probability that their Keplerian energy becomes negative. 

In order to investigate the influence of the capture boundary 
on the computed capture efficiency and to justify the 
plausibility of using three Hill Radii, we integrate all the 
particles within the circumsphere (seven Hill Radii) around P2 
for 6T = 2 yr. During the integration process for each particle, 
when the Keplerian energy is less than 0, the minimum distance 
of the particle from P2 is recorded. That is to say, particles with 
minimum distance from P2 less than three Hill Radii are 
precisely those that can be captured when the capture boundary 
is set to three Hill Radii. The distribution of particles’ minimum 
distances from P2 is displayed in Figure 11. According to the 
cumulative capture efficiency curve, it is evident that the 
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Figure 10. Linear fit for relation between log,(j) and log,)(p) based on uniform and non-uniform steady-state distributions. 
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Figure 11. The minimum distance from P2 for test particles with a negative Keplerian energy during the 2 yr integration. The histogram stands for the particle number. 
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capture efficiency increases as the capture boundary increases. 
Within one Hill Radius from P2, the number of particles 
decreases as the distance from P2 increases. This implies that 
particles have a higher likelihood of being found near P2 after 
entering the Hill Radius. Another peak in the distribution of the 
minimum distances of particles from P2 occurs around three to 
four Hill Radii. From an alternative perspective, when P2- 
centric distance is beyond one Hill Radius, particles with 
negative energy are more likely to be found in the range of 3-4 
Hill Radii. This may help explain why a three Hill Radii 
capture boundary serves as a choice for truncation of capture 
process. 


4.2. Orbit Eccentricity of P2 


The majority of our work is conducted under the framework 
of CRTBP model. The relative geometry between P1 and P2 
remains constant over time. The long-term gravitational 
influence of P2 can thus be investigated and the relationship 
between the capture efficiency and the mass parameter can be 
established only by altering the mass parameter. In the real 
context, planetary orbits are typically elliptic. In this subsec- 
tion, we extend our work to the elliptic restricted three-body 
problem (ERTBP) and alter the orbit eccentricity of P2 to 
examine whether it would influence the capture efficiency 
or not. 

The system used here is also the one with a mass parameter 
of the Sun—Earth CRTBP system. The orbit eccentricities of P2 
are listed in Table 2 as 0, 0.0167 (Earth-like orbit), 0.05 and 
0.1. The first thing is to determine the initial sampling region 
when the P2’s orbit eccentricity is non-zero. The same 
backward integration method in Section 2.3 is used both at 
periapsis and apoapsis of P2’s orbit. Our studies indicate that 
the size of the initial sampling region increases with P2’s orbit 
eccentricity. For the case of P2’s orbit eccentricity equals 0.1, 
the initial sampling region is identified as: a € [0.78, 1.28] au, 
e€[0, 0.17], i€ [0°, 2°5]. In order to compare the capture 
efficiency of different orbit eccentricity values at the same 
condition, we set the initial sampling region for all the tested 
orbit eccentricities as this one, and the bin resolution is: 
Aa = 0.01 au, Ae = 0.01, Ai=0°5. 

Basically, the same method as above has been used to 
compute the capture efficiency of each bin and the steady-state 
distribution. Considering that the relative geometry of ERTBP 
changes as P2 revolves around P1, an additional epoch is 
assigned to the particles during their generation. The epoch is 
uniformly distributed in the range of [0, 27] to average over the 
geometry between the particle, P2 and P1. Figure 12 displays 
the capture efficiency p(a, e, i) of each bin for different orbit 
eccentricity values. An obvious fact is that the center of the 
map remains close to the orbit eccentricity of P2. According to 
Table 2, the total capture efficiency decreases as the orbit 
eccentricity of P2 increases. This indicates that P2’s orbit 
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eccentricity has a negative influence on the total capture 
efficiency. It is noticeable that the capture efficiency for orbit 
eccentricity of P2 to be 0 in Figure 12 differs from the one 
computed in Section 3.2. This discrepancy is due to the fact 
that a larger initial sampling region is used here. It has been 
mentioned in Section | that the capture efficiency in our work 
serves as an indicator of P2’s ability to capture nearby particles 
within a specific region. The computed result is dependent on 
the choice of the initial sampling region and the alteration of 
the sampling region here is solely for the purpose of discussing 
the potential impact of P2’s orbit eccentricity on the capture 
efficiency. 


5. Conclusion 


In this paper, we study the temporary capture efficiency of 
massless particles around P2’s orbit in the framework of 
CRTBP. We focus on the gravitational influence of P2 on the 
capture efficiency. The final computed capture efficiency based 
on the non-uniform steady-state distribution which fully takes 
P2’s gravity into account is different from the one based on the 
uniform steady-state distribution. In the Sun—Earth CRTBP 
system (ip ~ 3.0035 x 1076), the capture efficiency based on 
the non-uniform distribution is computed as ~2.930 x 10-4. In 
this case, the initial sampling region which envelops P2’s orbit 
is given as a € [0.87, 1.15] au, e € [0, 0.06], i € [0°, 2°5]. The 
circumsphere size is set to 7 Hill Radii of P2. 

The way to compute the capture efficiency is the one used in 
Granvik et al. (2012). The relevant physical parameters during 
the process are also investigated and some reference values are 
recommended. A backward integration method is proposed to 
identify the initial sampling region. The circumsphere size 
around P2 is also altered to choose a proper value. With the 
relative increment in capture efficiency at specific circumsphere 
size less than 10%, 7 Hill Radii of P2 is recommended as the 
circumsphere size in the Sun—Earth CRTBP system. A 2 yr 
integration time is demonstrated to be long enough to compute 
the capture efficiency. Furthermore, we also find a power-law 
relation between the capture efficiency and mass parameter, 
which holds as 


p x 0.27 x pS 


when 3.0035 x 107% < u < 3.0034 x 1075. 

It is emphasized that since we do not consider realistic orbit 
element distributions of small bodies around the planet, the 
capture efficiency computed is not the real one of the planet, 
but it can serve as a metric for assessing the planet’s capability 
to capture nearby particles within a specific region. We believe 
that the planet’s gravity is a dominant factor which influences 
the steady-state distribution and the capture efficiency. The 
time-slicing method enables us to generate a steady-state 
distribution of nearby particles for a given planet. Moreover, 
although the algorithm is basically the same one taken by 
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Figure 12. Capture efficiency p(a, e, i) of each bin in the initial sampling region for different values of P2’s orbit eccentricity. The initial sampling region is a € [0.78, 1.28] au, 
e € [0, 0.17], i € [0°, 2°5] and the bin resolution remains: Aa = 0.01 au, Ae = 0.01, Ai = 0°5. The orbit eccentricities of P2 are: (1) Upper left: e = 0; (2) Upper right: 
e = 0.0167 (Earth-like orbit); (3) Lower left: e = 0.05; (4) Lower right: e = 0.1. The colorbar means the capture efficiency in the bin. 


Table 2 
Test Values of the ERTBP’s Orbit Eccentricity, Along with the Total Capture 
Efficiency 


0.0167 0.05 0.1 
p 5.928 x 1077 4.931 x 1075 2.142 x 107° 1.195 x 107" 


Orbit Eccentricity 0 


Granvik et al. (2012), the ways to choose parameters of the 
algorithm such as the size of the initial sampling region, 
circumsphere size and integration time can serve as good 
references if one wants to study the capture efficiency of some 
particular planet once a reliable orbit element distribution 
around it is available. The power-law relation between the 
capture efficiency and the mass parameter provides us with a 
quick reference of capture capability for any planets within the 
specific range (ug < 4<10pg). To be more realistic, we 
extend our work to the ERTBP. When the orbit eccentricity of 
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P2 increases, the computed capture efficiency decreases. This 
means that the orbit eccentricity of P2 has a negative influence 
on the capture efficiency. 

Our next step is to analyze the capture efficiency in the case 
where dissipative forces such as the Poynting—Robertson drag 
and the Yarkovsky effect are included. These thermal forces are 
significant for celestial bodies with small sizes. We try to find 
out how these thermal forces along with the gravity of P2 can 
affect the capture efficiency of particles in the vicinity. This 
may help us better understand the role of dissipative forces in 
the evolutionary history of small celestial objects around the 
planet, e.g., dust particles around Earth’s orbit. 
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